Electronic Pages

Die Homepage der Familie Beis

Analyse des Goertzel-Algorithmus

Erklärung, Zusatzinformationen, aber auch ein wenig Kritik

( English Version: Analysis of the Goertzel algorithm)

Inhalt

Ein paar Worte vorab

Hin und wieder fragt jemand in einem Forum, wie man eine spezielle Frequenz in einem digitalisierten Audiosignal erkennen kann. Als Antwort lese ich dann oft "Nimm doch einfach den Goertzel-Algorithmus" - oft genug mit dem Unterton "stell' dich doch nicht so dumm an!". Ich habe eine solche Anwendung bisher nur analog, aber noch nicht digital gehabt. Ich hätte einen digitalen IIR-Filter 2. Grades eingesetzt und so parametriert, dass er wie ein 2-poliger analoger Bandfilter mit ausgeprägter Resonanz arbeiten würde.

Aber mit dem Goertzel-Algorithmus schien es ja noch einfacher zu sein, und dumm wollte ich nicht sterben. Was ich im WWW bei meiner kurzen Suche dazu fand, war nicht sehr erhellend. Die Autoren wussten zwar offensichtlich, was sie meinten, aber ich, trotz meiner Kenntnisse von digitalen Filtern, eher nicht. Es sollte eine Sonderform der FFT sein, oder auch der DCT... Nun ja, irgendwann kam doch irgendwoher die Erleuchtung - und die Überraschung:

Der Goertzel-Algorithmus ist fast nichts anderes, als das, was ich mit einem digitalen IIR-Filter auch gemacht hätte. Und wegen der Unterschiede würde ich in der Praxis eher (m)einen IIR-Filter als den Goertzel-Algorithmus anwenden. Daher hätte die Antwort auf die Frage genauso gut (oder sogar noch besser) "Nimm doch einfach einen IIR-Filter" lauten können. Und wäre auch genauso hilfreich.

Kannte Gerald Goertzel schon die Theorie der digitalen Filter, als er 1958 diesen Algorithmus entwickelte?

Ich möchte hier nicht nur den Goertzel-Algorithmus erklären, sondern auch mit dem reinen IIR-Filter vergleichen, dazu das analoge Pendant. Und, was in den Foren offensichtlich nicht für erwähnenswert gehalten wird, auch die Parametrierung und andere zu beachtende Effekte der Filter erläutern. Letztendlich ist es viel mehr Text geworden, als ich vorgehabt hatte. Zu viel Text? Kann sein. Aber "Nimm doch einfach den Goertzel-Algorithmus" ist definitiv zu wenig Text.

Der ganz schnelle Einstieg

Dazu schreibe ich den Goertzel-Algorithmus als (nicht ganz vollständiges) Basic-Programm:

    Dim Z1, Z2 As Double        'Zwei interne Register mit größerer Wortbreite
    Dim X(BlockSize) as Integer 'Speicher für einen Block mit Samples
    SR = SampleRate             'Z. B. 48000 Hz
    FF = FilterFrequency        'Z. B. ein Rufton mit 1500 Hz
    N = BlockSize               'Anzahl der Samples, die zu einem Ergebnis führen
    Pi = 4 * Atn(1)             'Weil Basic Pi nicht kennt
    W = 2 * Pi * FF / SR        'Verhältnis FilterFrequenz/SampleRate * 2Pi
    CW = Cos(W)                 'Der Cosinus von W ist für höhere Frequenzen näher 1
    SW = Sin(W)                 'Der Sinus von W ist für höhere Frequenzen näher 0
    Do                          'Dauerschleife, je Durchlauf einen Block berechnen
        Z2 = 0                  'Die beiden internen Register werden zu Beginn
        Z1 = 0                  '    der Berechnung eines Blocks gelöscht
        For J = 1 To N          'Für alle N Samples eines Blocks:
            Z2 = Z1             '"Schieberegister"
            Z1 = Z0
            'X(J) ist der Wert dieses Samples:
            Z0 = X(J) + 2 * CW * Z1 - 1 * Z2 'Neuer Eingangswert des Schieberegisters
        Next J
        'Jetzt ist ein Block berechnet worden, das Ergebnis Y kann ausgegeben werden:
        Q = SW * Z1             'Der Realteil des Ergebnisses
        I = CW * Z1 - Z2        'Der Imaginärteil des Ergebnisses
        Y = Sqr(I ^ 2 + Q ^ 2)  'Y = Betrag (Länge des Vektors I,Q) ist das Ergebnis
    Loop

Ein analoger, äquivalenter Bandfilter

Ich schreibe hier über Bandfilter, aber eigentlich handelt es sich um 2-polige Tiefpassfilter mit ausgeprägter Resonanz am Bandende, also 2-polige Chebycheff-Tiefpassfilter mit sehr hoher Welligkeit. Ein solcher Tiefpassfilter kann als Analogschaltung aus Kapazität C, Induktivität L und Dämpfungswiderstand R bestehen und sieht z. B. wie folgt aus:

Wenn R klein genug ist, ergibt sich eine Resonanz bei der Frequenz, bei der die Scheinwiderstände X_L (von L) und X_C (von C) gleich sind. R bestimmt die Güte des Schwingkreises. Beispiel mit R = X_L/10 = X_C/10 = 10:

R = 10 Ω
L = 100 H (X_L = 100 Ω @ 159 mHz)
C = 10 mF (X_C = 100 Ω @ 159 mHz)

Niedrige Frequenzen werden 1:1 durchgelassen, hohe Frequenzen mit 12 dB/Oktave gedämpft, und bei der Resonanzfrequenz ergibt sich die 10-fache Amplitude, also die Güte Q = 10 (die Güte Q darf nicht mit dem Imaginärteil Q im obigen Programm verwechselt werden). Genau so kann sich auch ein normaler digitaler Filter verhalten, aber der im Goertzel-Algorithmus hat eine Eigenschaft, die mit analogen Filtern in der Realität nicht realisierbar ist: Er hat eine unendlich hohe Güte. Im analogen Pendant entspräche das R = 0 Ω und L und C sind ideale Bauelemente. Das lässt sich zumindest simulieren:

Bekanntlich pendelt in einem Schwingkreis Energie E immer zwischen zwei Energiespeichern, hier also zwischen der Kapazität mit der Energie E_C = 1 /2 · C · U² und der Induktivität mit der Energie E_L = 1 /2 · L · I². Die Gesamtenergie ist E = E_C + E_L. Sofern ein Eingangssignal präzise mit der Filterfrequenz übereinstimmt und dauerhaft anliegt, hat ein Filter mit unendlich hoher Güte die Eigenschaft, dass die Amplitude des Stroms und der Spannung im Laufe der Zeit linear und damit sein Energieinhalt quadratisch steigt. Genau das ist im Filter des Goertzel-Algorithmus der Fall. Dazu gleich mehr.

Auswertung des analogen Filterinhalts

Um nun mit so einem Filter mit extrem hoher Güte den spektralen Anteil seiner Frequenz im Signal zu messen, müsste der Filter für eine definierte Zeitspanne laufen, dann müsste als Messwert dessen Energieinhalt bestimmt werden, und für die nächste Messung müssten Spannung und Strom in L und C wieder auf 0 gesetzt werden. Das wäre für eine analoge Lösung zwar möglich, aber ziemlich unsinnig. Digital ist das einfacher, und im Goertzel-Algorithmus wird es genauso gemacht.

Der digitale Algorithmus

Schauen wir uns den digitalen Algorithmus als Blockschaltbild an:

Kurze Übersicht:

Der Algorithmus im linken grünen Bereich ist im Wesentlichen ein digitaler Filter 2. Grades in der Direktstruktur 2. Er entspricht der For-Next-Schleife im obigen Code sowie dem idealen analogen Filter (d. h., R = 0 und L und C ideal) und wird für jedes Eingangssample ausgeführt. Analog zur idealen analogen Schaltung steigt im Resonanzfall in seinen Registern das, was der Amplitude und der Energie entspricht, stetig. cw und sw stehen für den Cosinus bzw. Sinus des Verhältnisses von Abtastrate zu Filterfrequenz.

Zur Auswertung des Filterinhalts wird mit dem Algorithmus im rechten blauen Bereich nach einer definierten Anzahl von N Samples der bis dahin erreichte Energieinhalt ermittelt. Ausgegeben wird aber nicht die Energie I² + Q², sondern die Wurzel daraus, denn die Energie steigt mit dem Quadrat der Amplitude in den Registern, gewünscht wird aber ein Wert, der proportional zur Amplitude ist, also sqrt(I² + Q²).

Der digitale Bandfilter

Es gibt viele nicht uninteressante Zusammenhänge im digitalen Filter, zu dessen Anwendung man sie vielleicht nicht unbedingt verstehen muss, aber spätestens bei der Realisierung sollte man sich über einige Zusammenhänge bewusst sein, um nicht aus Unwissenheit "auf die Nase zu fallen".

Die unendliche Güte

Dass dieser Filter eine bestimmte Frequenz filtert und dass er eine unendliche Güte hat, lässt sich leicht plausibel machen:

• Nehmen wir den Sonderfall an, dass der Koeffizient a1 = 2·cw = 0 ist. In diesem Fall ist die zu Filterfrequenz = 1 / 4 der Abtastfrequenz.
• Außerdem nehmen wir an, dass sich die Eingangssamples alle 4 Takte wiederholen, z. B.: 1  0 ‑1  0 oder auch 1 1 0 0, also ebenfalls 1 / 4 der Abtastfrequenz:
• Dabei wird das rückgeführte Signal nach jeweils 2 Takten in derselben Phase bzw. Polarität zum Eingangssignal addiert, wodurch die Registerinhalte in jedem Zyklus um 1 erhöht werden.
• Sofern sie nicht überlaufen, ist der erreichte Wert in den Registern proportional zu der Zeit, die dieser Prozess läuft. Wenn es einen idealen analogen Filter gäbe, wäre das bei dem genauso.

Wegen der begrenzten Anzahl Samples bedeutet eine unendliche Güte aber nicht, dass nur Spektralanteile mit genau der theoretischen Filterfrequenz detektiert werden. Die begrenzte Anzahl Samples wirkt sich in der Praxis genauso aus wie ein Filter mit endlicher Güte. Das bedeutet, dass auch Spektralanteile nahe der Filterfrequenz in die Messung, lediglich mit entsprechend reduzierter Amplitude, einfließen.

Filterfrequenz

Sofern, wie beim Goertzel-Algorithmus, der Koeffizient a2 = ‑1 ist, bestimmt der Koeffizient a1 = 2·cw = 2 · cos(2 · Pi · FF / SR) die Resonanz- bzw. Filterfrequenz. Er ist < 2, und je näher er an 2 liegt, desto niedriger ist die Filterfrequenz.

Beispiel: Es soll eine Frequenz von 1500 Hz bei einer Abtastfrequenz von 48000 Hz berechnet werden. Es ergibt sich: a1 = 2·cw = 2 · cos(2 · Pi · 1500 / 48000) = 1,96157056080646

Nebenbei: Nicht nur die Filterfrequenz wird durch den Koeffizienten a1 bestimmt, sondern auch die DC-Verstärkung G wird verändert: G = 0.5 / (1 ‑ cw). Im vorherigen Beispiel wäre G = 0.5 / (1 ‑ cw) = 26.0217172299543, wodurch alle Spektralanteile unterhalb der Resonanzfrequenz um mindestens den Faktor 26 verstärkt in den Registerwerten erscheinen würden. In der Praxis ist das aber eher unerheblich, denn je niedriger die Resonanzfrequenz ist, desto höher ist zwar dieser Einfluss, aber gleichzeitig ist der Energieinhalt im Spektrum unterhalb der Resonanzfrequenz auch geringer.

Experimente mit AnaDigFilt

Ich habe vor langer Zeit das Programm AnaDigFilt geschrieben, mit dem digitale Filter aus analogen Vorlagen dimensioniert werden können. Es lassen sich aber auch Frequenzgänge digitaler Filter anhand gegebener Koeffizienten darstellen. Das war beim Schreiben dieses Artikels sehr hilfreich. Es ist (leider) eine EXE-Datei für Windows, die man auf der Seite Converting Analog into Digital (IIR) Filters herunterladen kann. Hinweis: Um den Frequenzgang des Beispiels darzustellen, muss eingetragen werden:

Im Kopf die analogen Koeffizienten:
    Fc = 1500 Hz (nur, damit eine Linie bei 1,5 kHz erscheint)
    Filter Order n = 2 (weil es ein Filter 2. Grades ist)
Die digitalen Koeffizienten werden nicht aus den analogen berechnet, sondern manuell eingesetzt:
    entweder die oberen aus der oberen Gleichung dieses Screenshots oder die unteren aus der unteren Gleichung.
Dann muss noch ein geeigneter Diagrammausschnitt eingestellt werden.

Die Resonanzfrequenz von 1,5 kHz und der DC-Verstärkungsfaktor von ca. 26 = ca. 28 dB ist nachvollziehbar.

Registerbreite bei der Berechnung

Bei jedem digitalen Filter gibt es einen Zusammenhang, von dem ich fürchte, dass er leicht übersehen wird, wenn man den Hintergrund nicht kennt. Das kann zu schlimmen Konsequenzen führen: Der Eingangswert, der eine bestimmten Wortbreite bzw. Anzahl Bits bzw. Auflösung hat, sollte bei der folgenden Berechnung nicht unnötig reduziert werden. Bei Filtern kann das aber schnell passieren:

• Je niedriger die Filterfrequenz ist, desto mehr Samples werden benötigt, um die gleiche Filtergüte Q zu bekommen.
• Je mehr Samples verarbeitet bzw. akkumuliert werden, also je größer N ist, desto mehr Stellen werden bei der Berechnung benötigt.
• Dabei entstehen oft sehr kleine Differenzen aus zwei sehr großen Werten, deswegen dürfen die vielen Stellen nicht reduziert bzw. gerundet werden.

Berücksichtigt man das nicht, kann der berechnete Wert bis zur völligen Unbrauchbarkeit verfälscht werden.

Beispiel: Bei einem Filter mit einer angenommenen Filterfrequenz von 50 Hz und der Abtastfrequenz von 48000 Hz sowie einer gewünschten Güte Q = 100 (das entspricht einer Bandbreite von 0,5 Hz) werden 100·48000/50 = 96000 Samples (2 Sekunden) benötigt. Also können die Werte in den Registern fast das 100000-fache des Eingangssignals annehmen, sodass es entweder mehr Vorkommastellen gibt, oder das Eingangssignal wird mit Koeffizienten k = 1/96000 multipliziert (sinnvollerweise durch 2^17 = 131072 geteilt, damit wird nur eine Shift-Operation mehr benötigt), was entsprechend viele Nachkommastellen erfordert. In jedem Fall sind in diesem Beispiel zusätzliche 17 Bit bei der Berechnung erforderlich.

Lägen in diesem Fall Samples am Eingang als 16-Bit-Werte (Integer) vor und würde die folgende Rechnung auch nur mit 16 Bit ausgeführt, würden selbst bei sehr kleinen Eingangswerten keine sinnvollen Ausgangswerte mehr entstehen.

Anderes Beispiel: Soll die Filterfrequenz wieder 50 Hz sein, ergibt sich der Koeffizient a1 = 2·cw = 1,99995716332826. Im Nenner der Gleichung H(z) steht dann 1 ‑ 1,99995716332826 + 1 = 0,0000428366717399875. Den Wert von a1 jetzt zu 1,99996 aufzurunden - das ist ja fast dasselbe - ergäbe eine Resonanzfrequenz von ca. 48 Hz. Hier sogar mit 32-Bit float zu rechnen, wäre "ganz dünnes Eis".

Konsequenz:

Weil die Werte in den Registern maximal das N-fache des Eingangssignals ergeben können, werden sinnvollerweise die Eingangssamples oder alternativ später auch die Ausgangssignale durch N geteilt, um den gleichen Wertebereich wie den des Eingangssignals zu erhalten. Außerdem sollte mit einer entsprechend hohen Auflösung gerechnet werden, was insbesondere bei niedrigen Grenz- bzw. Resonanzfrequenzen sehr viele Stellen umfassen kann.

Auswertung des digitalen Filterinhalts

Die Energie bzw. das Äquivalent der Energie, das in einem digitalen Filter vorliegt, lässt sich berechnen aus

• dem momentanen Inhalt eines der Register als Realteil (vergleichbar mit der Ladung C · U im Kondensator) und
• der Änderungsgeschwindigkeit des Registerinhalts als Imaginärteil (vergleichbar mit L · I, wofür es keinen Begriff wie "Ladung" beim Kondensator gibt).

Das Maximum, das für Y über eine Messperiode entstehen kann, entspricht dem spektralen Anteil, der im Signal mit der Resonanzfrequenz enthalten ist, akkumuliert über die Anzahl der Samples. Im Idealfall wird für ein Eingangssignal mit der Amplitude = 1 nach N Samples Y = 0,5 · N entstehen.

Darstellung im Zeitbereich

Grafisch lässt sich das Verhalten des Goertzel-Algorithmus veranschaulichen. Bei sinusförmigen Eingangssignalen mit der Amplitude 1 werden für 3 verschiedene Frequenzen über 1200 Samples die Werte Z0 und die bis dahin erreichten Messwerte Y dargestellt:

• Blau: Filter- und Signalfrequenz 1,5 kHz: Die Amplitude im Filter wächst kontinuierlich. Ebenso wächst die Energie im Filter und die mit der oberen rote Linie dargestellte Quadratwurzel Y daraus, bis der Messwert Y nach 1000 Samples den Wert 500 erreicht.
• Gelb: Bei einer anderen, hier z. B. ungefähr der doppelten Filter- und Signalfrequenz 3,2 kHz, wächst die Amplitude im Filter nur ungefähr halb so schnell (siehe Registerbreite), aber der Messwert Y wächst genauso schnell (ebenfalls die obere rote Linie, die sich mit der vorherigen deckt).
• Grün: Die Signalfrequenz liegt 5% über der Filterfrequenz. Dabei entsteht nur zu Anfang eine etwas "pulsierende" Amplitude, die Energie bzw. der Messwert Y schwankt und bleibt dauerhaft klein (untere rote Linie).

Nach N = 1000 Samples werden die Register gelöscht, und der Messvorgang wird neu gestartet.

Weshalb ich einen fast identischen, einfachen IIR-Filter bevorzuge

Beim Goertzel-Algorithmus werden Daten, wie bei einer DCT oder FFT, immer blockweise berechnet und am Ende des Blocks wird ein Ergebnis geliefert. Dadurch werden z. B. für kurze Bursts, je nachdem wie ihre jeweilige zeitliche Lage zu der der Blöcke ist, unterschiedliche Messwerte ausgegeben. Bei einer DCT oder FFT ist das genauso.

Die Messwertausgabe kann nur in regelmäßigen Zeitabständen springen, ähnlich, wie es bei digitalen Multimetern der Fall ist, und das gefällt mir nicht. Aber im Gegensatz zu digitalen Multimetern bieten sie nicht den Vorteil höherer Auflösung. Außerdem wird die Messwertausgabe ständig mehr oder weniger springen, wenn eine Frequenz neben der exakten Resonanzfrequenz vorliegt.

Deswegen ist mir das Verhalten eines analogen Filters, wie auch das Verhalten eines "normalen" digitalen IIR-Filters mit endlicher Güte, sympathischer: Mit jedem Sample entsteht in den Registern ein Wert, den ich jederzeit und nicht nur in größeren Zeitabständen, z. B. auf die Überschreitung eines Schwellwerts oder zur kontinuierlichen, quasi-analogen Anzeige, verwenden kann. Nachteil: Je nach Anforderung bzw. "Perfektion" der kontinuierlichen Ausgabe ist mehr oder weniger zusätzlicher Rechenaufwand erforderlich:

Abgesehen von dem Faktor a2 im Signalweg Z2, der für eine endliche Güte zwischen 0 und ‑1 liegt, wäre der Filter identisch zum grünen Bereich des Goertzel-Algorithmus. Zur Auswertung auf die Überschreitung einer Schwelle braucht lediglich ein Registerwert auf die Überschreitung eines Schwellwerts geprüft zu werden. Für eine kontinuierliche, quasi-analoge Anzeige, wäre der rechte blaue Teil des Goertzel-Algorithmus die mathematisch perfekte Lösung, es ginge aber auch einfacher mit einem Gleichrichter (Absolutwert) und einem nachgeschalteten einfachen Tiefpass.

IIR-Bandfilter mit Q = 30 im Zeitbereich

Zum Vergleich zeige ich das Verhalten eines "normalen" digitalen IIR-Bandfilters mit endlicher Güte Q = ca. 30. Die Parameter a1 = 1,8145, a2 = ‑0,987 und k = 0,173 habe ich nicht exakt berechnet, sondern mit AnaDigFilt experimentell bestimmt. Damit die Skalierung zur vorherigen passt, ist hier die Amplitude der Anregung allerdings = 5,4. Nach ca. 11 Schwingungen (3,7 ms) erreichen die Amplituden in den Registern ca. 70% der endgültigen Amplitude von ca. 1000.

Auch hier schwingt sich der Filter bei etwas abseits liegender Signalfrequenz erst etwas stärker ein, bleibt aber auf die Dauer kleiner. Die Welligkeit der Ausgabe Y (rote Linie) wird durch die Dämpfung des Filters verursacht. Im Gegensatz zum idealen Filter verliert dieser verlustbehaftete Filter zwischen den Halbwellen etwas Energie und gewinnt sie während jeder Halbwelle wieder.

Tipp: a2 = ‑0,987 muss nicht unbedingt präzise eingehalten werden. a2 ist entscheidend für die Dämpung bzw. Güte des Filters, die muss oft nicht sehr genau stimmen. Wenn die Ressourcen der Hardware knapp sind, kann man sparen: Mit a2 = ‑1 + 1/128 + 1/256, also mit nur zwei zusätzlichen Shift-/Add-Operationen, wird in diesem Beispiel schon ein in der Praxis oft ausreichender Wert für a2 von 0,9882813 erreicht. Mit einer weiteren Shift/Add-Operation a2 = ‑1 + 1/128 + 1/256 + 1/1024 = 0,9873047 ist praktisch kein Unterschied in der Verhaltensweise des Filters mehr zu erkennen.

DTMF-Dekodierung mit dem Goertzel-Algorithmus

Das DTMF-Wählverfahren brauche ich hier sicherlich nicht zu erklären. In der Praxis dürfte der Goertzel-Algorithmus das am häufigsten angewendete Verfahren zu Erkennung der DTMF-Töne sein. Allerdings besteht ein DTMF-Dekoder nicht nur aus 8 Dekodern für die 8 Töne. Sie müssen geprüft werden: Signaldauer und Signalstärke sowie der Twist bzw. der Reverse Twist müssen im Toleranzbereich liegen, damit ein Tonpaar als gültig bzw. ungültig erkannt werden kann. Schlussendlich muss ein gültiges Tonpaar einen von 16 Ausgängen aktivieren.

Soll die volle Signalstärke in einer Messperiode erkannt werden, dürfte gemäß den Spezifikationen die Messperiode höchstens 20 ms lang sein. Andererseits sollte es möglich sein, auch die Summe der Signalstärken von zwei aufeinander folgenden, 40 ms langen Messperioden auszuwerten. Das würde zu geringerer Filterbreite, aber auch zu langsameren Reaktionen führen und auch längere Pausen zwischen zwei gleichen Tönen erfordern. Eigentlich sollen auch bei Pausen im Signal von mehr als 10 ms schon zwei getrennte Signaltöne erkannt werden. Das ginge nur mit wesentlich geringeren Messzeiten und entsprechend größeren Filter-Bandbreiten.

Aber alles das will ich hier gar nicht genauer diskutieren oder gar eine komplette Lösung vorstellen. Ich will damit nur andeuten, dass "den Goertzel-Algorithmus anwenden" nur ein Teil der Lösung für einen DTMF-Dekoder ist. Dieser Artikel soll lediglich helfen, die Parameter für die Filter eines DTMF-Dekoders zu bestimmen. Hier gibt es einen Vorteil bei einem Filter mit unendlicher Güte: Nur ein einziger, von Abtast- und Filterfrequenz abhängiger Filter-Parameter ist für die Filterfrequenz zuständig:

a1 = 2 · cos(2 · Pi · FF / SR)

Dazu zur Auswertung:

cw = a1 / 2
sw = sin(2 · Pi · FF / SR)

Für einen IIR-Filter mit begrenzter Güte müsste auch der Parameter a2 berechnet werden. Das ist komplexer, denn beide Parameter, a1 und a2, sind abhängig von der Güte und der Abtast- und der Filterfrequenz. Ich habe aber keine Formeln, um das zu berechnen. Dennoch, ich würde für einen DTMF-Dekoder wahrscheinlich diesen Weg gehen.